Objectives

Introduction

So we’ve made it this far into the quarter and it’s time for the last lab. There is not final exam but this lab will consist of many of the past lessons from the previous 9. By now you should have a repository of code and lidar processing examples in the form of your past assignments. Most of this final lab comes from those past assignments and reuses their code. So there will be more limited examples here meaning less copy/paste.

The data in this lab also consists of past lab examples except for some extra TLS data from Pack Forest and the Hall of Mosses in the Hoh Rainforest: 2019-08-18 Pack Forest TMLS.laz& 2019-08-22 Hall of Mosses TMLS.laz. You will need to download extra data from the lidar portal in Part 1.

Part 1: ALS & TLS comparison processing steps with required figures

You will have 2 plots for this part centered at:

Download the ALS point clouds at those two locations from the Washington DNR lidar portal website.

PFTLS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/2019-08-18 Pack Forest TMLS.laz')
HMTLS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/2019-08-22 Hall of Mosses TMLS.laz')
PFALS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/WA_031_rp.laz')
HMALS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/q47123g8201.laz')
NEWCRS <- crs("EPSG:2927")
crs(PFTLS) <- NEWCRS
crs(HMTLS) <- NEWCRS
crs(PFALS) <- NEWCRS
crs(HMALS) <- NEWCRS
?clip_circle
PFTLSc <- clip_circle(PFTLS,1185950,555180,50)
HMTLSc <- clip_circle(HMTLS,797785,941189,50)
PFALSc <- clip_circle(PFALS,1185950,555180,50)
HMALSc <- clip_circle(HMALS,797785,941189,50)

Future steps for the TLS will only use these clips

writeLAS the four lidar clips and bring them into CloudCompare to produce a figure.

?writeLAS
writeLAS(PFTLSc)
writeLAS(HMTLSc)
writeLAS(PFALSc)
writeLAS(HMALSc)

REQUIRED FIGURES: An image of the ALS and TLS plotted together for both the Pack Forest (PF) and Hall of Mosses (HM) sites (two screen shots). Example given has the TLS cloud in gray scale and the ALS cloud in multicolor for the Hall of Mosses forest plot.

Now back in R:

mycsf <- csf(FALSE, 1, 1, rigidness = 3)
#OR
mypmf <- pmf(ws = 0.5, th = 0.1)

PFTLSg <- classify_ground(PFTLSc, mycsf)

HMTLSg <- classify_ground(HMTLSc, mycsf)

plot(PFTLSg, color = "Classification")
PFDTM1 <- rasterize_terrain(PFTLSg, res = 1, algorithm = knnidw())

HMDTM1 <-grid_terrain(HMTLSg, res = 1, algorithm = knnidw())
## Warning: There were 1 degenerated ground points. Some X Y Z coordinates were
## repeated. They were removed.
## Warning: There were 172 degenerated ground points. Some X Y coordinates were
## repeated but with different Z coordinates. min Z were retained.
plot_dtm3d(PFDTM1)
plot_dtm3d(HMDTM1)

REQUIRED FIGURES: Two screen shots of your two produced dtms. Fully captioned with the technique used to make them.

PFTLScN <- normalize_height(PFTLSc,PFDTM1)
HMTLScN <- normalize_height(HMTLSc,HMDTM1)
PFALScN <- normalize_height(PFALSc,PFDTM1)
HMALScN <- normalize_height(HMALSc,HMDTM1)
?rasterize_canopy
PFTLS_DSM <- rasterize_canopy(PFTLScN,res = 1.5, p2r(2))
PFALS_DSM <- rasterize_canopy(PFALScN,res = 1.5, p2r(2))  
HMTLS_DSM <- rasterize_canopy(HMTLScN,res = 1.5, p2r(2))
HMALS_DSM <- rasterize_canopy(HMALScN,res = 1.5, p2r(2)) 
plot_dtm3d(PFTLS_DSM)
plot(PFALScN)
plot(HMTLScN)
plot(HMALScN)
?rumple_index
rumple_index(PFTLS_DSM)
## [1] 4.929346
rumple_index(PFALS_DSM)
## [1] 5.576009
rumple_index(HMTLS_DSM)
## [1] 3.509282
rumple_index(HMALS_DSM)
## [1] 5.997609

REQUIRED FIGURES: Four screen shots. One of each of your produced DSMs with the rumple index included in the caption along with the algorithm used to produce the DSM.

- Use focal statistics to smooth your DSM. You can do this step before generating the screenshots and rumple index but the smoothing may cause an error with the rumple index.

?focal
## Help on topic 'focal' was found in the following packages:
## 
##   Package               Library
##   raster                /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
##   terra                 /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
## 
## 
## Using the first match ...
PFTLS_DSMs <- focal(PFTLS_DSM, w = matrix(1,3,3), fun = mean)
PFALS_DSMs <- focal(PFALS_DSM, w = matrix(1,3,3), fun = mean)       
HMTLS_DSMs <- focal(HMTLS_DSM, w = matrix(1,3,3), fun = mean)
HMALS_DSMs <- focal(HMALS_DSM, w = matrix(1,3,3), fun = mean) 
plot_dtm3d(PFTLS_DSMs)
plot_dtm3d(PFALS_DSMs)
plot_dtm3d(HMTLS_DSMs)
plot_dtm3d(HMALS_DSMs)
PFTLStrees <- segment_trees(PFTLScN, lidR::watershed(PFTLS_DSMs, th = 10))
plot(PFTLStrees,color = "treeID", pal = pastel.colors)
PFALStrees <- segment_trees(PFALScN, lidR::watershed(PFALS_DSMs, th = 10))
plot(PFALStrees,color = "treeID", pal = pastel.colors)
HMTLStrees <- segment_trees(HMTLScN, lidR::watershed(HMTLS_DSMs, th = 10))
plot(HMTLStrees,color = "treeID", pal = pastel.colors)
HMALStrees <- segment_trees(HMALScN, lidR::watershed(HMALS_DSMs, th = 10))
plot(HMALStrees,color = "treeID", pal = pastel.colors)
#OR
PFTLSttops <- locate_trees(PFTLS_DSMs, lmf(ws = 3.28, hmin = 6.56))
algo1 <- dalponte2016(PFTLS_DSMs, PFTLSttops, th_tree = 0.1, th_seed= 0.1, max_cr = 10, th_cr = 0.1)
algo2 <- silva2016(PFTLS_DSMs, PFTLSttops, max_cr_factor = 0.6, exclusion = 0.3)


PFTLStrees1 <- segment_trees(PFTLScN, algo1)
plot(PFTLStrees1,color = "treeID", pal = pastel.colors)
PFTLStrees2 <- segment_trees(PFTLScN, algo2)
plot(PFTLStrees2,color = "treeID", pal = pastel.colors)

PFALStrees2 <- segment_trees(PFALScN, lidR::watershed(PFALS_DSMs, th = 10))
plot(PFALStrees,color = "treeID", pal = pastel.colors)
HMTLStrees <- segment_trees(HMTLScN, lidR::watershed(HMTLS_DSMs, th = 10))
plot(HMTLStrees,color = "treeID", pal = pastel.colors)
HMALStrees <- segment_trees(HMALScN, lidR::watershed(HMALS_DSMs, th = 10))
plot(HMALStrees,color = "treeID", pal = pastel.colors)

REQUIRED FIGURES: Four screen shots: for the ALS & TLS in each plot of the lidar data colored by the tree segmentation for each of the normalized lidar clips. Fully captioned with the algorithms used to create the tree segmentation.

Calculate cloud metrics on the normalized point cloud for the ALS and TLS above 2m height (6.56ft)

?cloud_metrics
PFTLS_cm <- cloud_metrics((filter_poi(PFTLScN,Z > 6.56)), ~stdmetrics_z(Z))
PFTLS_cm
## $zmax
## [1] 140.6
## 
## $zmean
## [1] 25.78389
## 
## $zsd
## [1] 28.6307
## 
## $zskew
## [1] 1.809775
## 
## $zkurt
## [1] 4.91211
## 
## $zentropy
## [1] 0.7387645
## 
## $pzabovezmean
## [1] 22.80502
## 
## $pzabove2
## [1] 100
## 
## $zq5
## [1] 7
## 
## $zq10
## [1] 7.5
## 
## $zq15
## [1] 8
## 
## $zq20
## [1] 8.5
## 
## $zq25
## [1] 9
## 
## $zq30
## [1] 9.7
## 
## $zq35
## [1] 10.4
## 
## $zq40
## [1] 11.1
## 
## $zq45
## [1] 12
## 
## $zq50
## [1] 12.9
## 
## $zq55
## [1] 13.9
## 
## $zq60
## [1] 15.3
## 
## $zq65
## [1] 16.9
## 
## $zq70
## [1] 19
## 
## $zq75
## [1] 23
## 
## $zq80
## [1] 32.2
## 
## $zq85
## [1] 61.8
## 
## $zq90
## [1] 80.6
## 
## $zq95
## [1] 97.4
## 
## $zpcum1
## [1] 55.57681
## 
## $zpcum2
## [1] 78.54518
## 
## $zpcum3
## [1] 82.06199
## 
## $zpcum4
## [1] 83.66011
## 
## $zpcum5
## [1] 86.83228
## 
## $zpcum6
## [1] 91.50673
## 
## $zpcum7
## [1] 95.48767
## 
## $zpcum8
## [1] 99.00548
## 
## $zpcum9
## [1] 99.97011
PFALS_cm <- cloud_metrics((filter_poi(PFALScN,Z > 6.56)), ~stdmetrics_z(Z))
PFALS_cm
## $zmax
## [1] 148.57
## 
## $zmean
## [1] 90.02078
## 
## $zsd
## [1] 42.76528
## 
## $zskew
## [1] -0.9650657
## 
## $zkurt
## [1] 2.296773
## 
## $zentropy
## [1] 0.8806774
## 
## $pzabovezmean
## [1] 69.47884
## 
## $pzabove2
## [1] 100
## 
## $zq5
## [1] 10.04
## 
## $zq10
## [1] 15.74
## 
## $zq15
## [1] 19.4475
## 
## $zq20
## [1] 25.13
## 
## $zq25
## [1] 73.0275
## 
## $zq30
## [1] 89.14
## 
## $zq35
## [1] 96.9475
## 
## $zq40
## [1] 102.82
## 
## $zq45
## [1] 106.3125
## 
## $zq50
## [1] 109.495
## 
## $zq55
## [1] 111.9375
## 
## $zq60
## [1] 114.04
## 
## $zq65
## [1] 116.15
## 
## $zq70
## [1] 118.59
## 
## $zq75
## [1] 120.63
## 
## $zq80
## [1] 122.53
## 
## $zq85
## [1] 125.02
## 
## $zq90
## [1] 127.975
## 
## $zq95
## [1] 131.5575
## 
## $zpcum1
## [1] 8.970482
## 
## $zpcum2
## [1] 22.23182
## 
## $zpcum3
## [1] 23.87329
## 
## $zpcum4
## [1] 24.13247
## 
## $zpcum5
## [1] 25.19798
## 
## $zpcum6
## [1] 30.0216
## 
## $zpcum7
## [1] 41.79986
## 
## $zpcum8
## [1] 70.59755
## 
## $zpcum9
## [1] 96.74586
HMTLS_cm <- cloud_metrics((filter_poi(HMTLScN,Z > 6.56)), ~stdmetrics_z(Z))
HMTLS_cm
## $zmax
## [1] 110.178
## 
## $zmean
## [1] 28.15648
## 
## $zsd
## [1] 20.19887
## 
## $zskew
## [1] 1.026229
## 
## $zkurt
## [1] 3.162441
## 
## $zentropy
## [1] 0.8600562
## 
## $pzabovezmean
## [1] 37.9027
## 
## $pzabove2
## [1] 100
## 
## $zq5
## [1] 7.292
## 
## $zq10
## [1] 8.09
## 
## $zq15
## [1] 9.077
## 
## $zq20
## [1] 10.4
## 
## $zq25
## [1] 11.863
## 
## $zq30
## [1] 13.384
## 
## $zq35
## [1] 14.856
## 
## $zq40
## [1] 16.642
## 
## $zq45
## [1] 18.863
## 
## $zq50
## [1] 21.21
## 
## $zq55
## [1] 23.4359
## 
## $zq60
## [1] 26.506
## 
## $zq65
## [1] 30.843
## 
## $zq70
## [1] 36.105
## 
## $zq75
## [1] 40.834
## 
## $zq80
## [1] 45.895
## 
## $zq85
## [1] 52.703
## 
## $zq90
## [1] 60.056
## 
## $zq95
## [1] 69.008
## 
## $zpcum1
## [1] 22.20074
## 
## $zpcum2
## [1] 51.81598
## 
## $zpcum3
## [1] 66.94687
## 
## $zpcum4
## [1] 78.20104
## 
## $zpcum5
## [1] 86.71497
## 
## $zpcum6
## [1] 93.61396
## 
## $zpcum7
## [1] 97.68797
## 
## $zpcum8
## [1] 99.47636
## 
## $zpcum9
## [1] 99.97858
HMALS_cm <- cloud_metrics((filter_poi(HMALScN,Z > 6.56)), ~stdmetrics_z(Z))
HMALS_cm
## $zmax
## [1] 111.88
## 
## $zmean
## [1] 68.30136
## 
## $zsd
## [1] 23.11787
## 
## $zskew
## [1] -0.7348024
## 
## $zkurt
## [1] 2.969791
## 
## $zentropy
## [1] 0.9386929
## 
## $pzabovezmean
## [1] 56.63975
## 
## $pzabove2
## [1] 100
## 
## $zq5
## [1] 20.18
## 
## $zq10
## [1] 34.552
## 
## $zq15
## [1] 43.418
## 
## $zq20
## [1] 49.264
## 
## $zq25
## [1] 54.43
## 
## $zq30
## [1] 59.81
## 
## $zq35
## [1] 63.262
## 
## $zq40
## [1] 66.23
## 
## $zq45
## [1] 69.464
## 
## $zq50
## [1] 72.74
## 
## $zq55
## [1] 75.05
## 
## $zq60
## [1] 77.552
## 
## $zq65
## [1] 80.428
## 
## $zq70
## [1] 83.04
## 
## $zq75
## [1] 85.66
## 
## $zq80
## [1] 88.512
## 
## $zq85
## [1] 91.902
## 
## $zq90
## [1] 95.058
## 
## $zq95
## [1] 98.868
## 
## $zpcum1
## [1] 2.131886
## 
## $zpcum2
## [1] 5.826271
## 
## $zpcum3
## [1] 9.666314
## 
## $zpcum4
## [1] 16.19439
## 
## $zpcum5
## [1] 26.20498
## 
## $zpcum6
## [1] 41.56515
## 
## $zpcum7
## [1] 61.28178
## 
## $zpcum8
## [1] 81.43538
## 
## $zpcum9
## [1] 96.875
PFTLS_cma <- cloud_metrics(PFTLScN, ~sum(Z>6.56)/sum(Z>-1))
PFTLS_cma
## [1] 0.5067569
PFALS_cma <- cloud_metrics(PFALScN, ~sum(Z>6.56)/sum(Z>-1))
PFALS_cma
## [1] 0.963117
HMTLS_cma <- cloud_metrics(HMTLScN, ~sum(Z>6.56)/sum(Z>-1))
HMTLS_cma
## [1] 0.5147704
HMALS_cma <- cloud_metrics(HMALScN, ~sum(Z>6.56)/sum(Z>-1))
HMALS_cma
## [1] 0.6470487
name Zmax Zmean Zsd Zq25 Zq95 CanopyCover
Pack Forest ALS 148.570 90.02078 42.76528 73.0275 131.5575 0.9631170
Pack Forest TLS 140.600 25.78389 28.63070 9.0000 97.4000 0.5067569
ONP ALS 111.880 68.30136 23.11787 54.4300 98.8680 0.6470487
ONP TLS 110.178 28.15648 20.19887 11.8630 69.0080 0.5147704

That is all for the processing of the ALS/TLS plot

Part 2: ALS & GEDI comparison processing steps with required figures for the report

This part is dealing with some more recent coding so I’ll provide it here. The questions will refer to the outputs

You will be comparing ALS cloudmetrics and Canopy Height Models (CHMs) to GEDI for the Hall of Mosses area of the Olympic National Park

For the purpose of this lab, we are going to assume that there is no error in the GEDI pulse locations. In reality there is some error and that can account for oddities in the data.

To complete this portion, of the lab you will need to:

Download The DTM and DSM from the Olympic Peninsula

Read them into R and assign the ALS the correct crs

ONP_DTM <- rast('Lab 10/Lab_10_2023/LAB10_Data_2023/hoh_2013_dtm_4.tif')
ONP_DSM <- rast ('Lab 10/Lab_10_2023/LAB10_Data_2023/hoh_2013_dsm_4.tif')

Subtract the DTM from the DSM for an easy Canopy Height Model

#ONPCHM <- ONPDSM-ONPDTM
#plot(ONPCHM)
ONPCHM <- rast('Lab 10/Lab_10_2023/LAB10_Data_2023/out/ONPCHM.tif')
#writeRaster(ONPCHM, ".....................ONPCHM.tif")

Now read in GEDI level 2a data from the data folder for this lab

?getLevel2AM
gedilevel2a<-readLevel2A(level2Apath = file.path ("Lab 10/Lab_10_2023/LAB10_Data_2023/HM_level2a_clip.h5"))

Compute a series of statistics of GEDI Level2A RH100 metric

#Get GEDI Elevation and Height Metrics (GEDI Level2A)
?getLevel2AM
level2AM<-getLevel2AM(gedilevel2a)
head(level2AM[,c("beam","shot_number","elev_highestreturn","elev_lowestmode","rh100")])

Convert shot_number as “integer64” to “character” and the make the GEDI Level2AM points into a spatial object with sf

level2AM$shot_number<-as.character(level2AM$shot_number)
level2AM_sf<-st_as_sf(level2AM, coords = c("lon_lowestmode", "lat_lowestmode"), crs = 4326)

Define a function for computing GEDI stats

mySetOfMetrics = function(x)
{
  metrics = list(
    min =min(x), # Min of x
    max = max(x), # Max of x
    mean = mean(x), # Mean of x
    sd = sd(x)# Sd of x
  )
  return(metrics)
}

Create a RH100 raster with 0.005 res in degrees of lat / long.

?gridStatsLevel2AM
rh100metrics<-gridStatsLevel2AM(level2AM = level2AM, func=mySetOfMetrics(rh100), res=0.005)
crs(rh100metrics) <- "EPSG:4326" #same as GEDI
rh100metrics_prj <- terra::project(rast(rh100metrics), "EPSG:2927") #transform to WA state plane south

plot(ext(770000, 850000, 900000, 970000), type = "lines", buffer = TRUE)
plot(rh100metrics_prj$max, axes = F, add = TRUE, legend =T)

plot(ext(770000, 850000, 900000, 970000), type = "lines", buffer = TRUE)
plot(ONPCHM, axes = F, add = TRUE, legend =T)

writeRaster(rh100metrics_prj, "Lab 10/Lab_10_2023/LAB10_Data_2023/out/rh100metrics_prj.tif")

REQUIRED FIGURES: Create 2 figures: a map from the ONPCHM and a map of the same area using the GEDI RH100 max data. Write out the ONPCHM and GEDI RH100 to open in ArcGIS Pro and put them next to each other for comparison and qualitatively assess the spatial patterns. In your caption explain if they generally agree or disagree with each other

Extract ONPCHM elevation values with the new projected GEDI Level2A points

ONPCHM <- rast("Lab 10/Lab_10_2023/LAB10_Data_2023/out/ONPCHM.tif")
level2AGeo_sf_prj <- st_transform(level2AM_sf, "EPSG:2927") #transform to WA state plane south
ONP_CHM_GEDI<- (terra::extract(ONPCHM, level2AGeo_sf_prj, bind = TRUE))
ONP_CHM_GEDI_FILTER <- ONP_CHM_GEDI[-c(is.na(ONP_CHM_GEDI$hoh_2013_dsm_4)),] #removing NA values 
ONP_CHM_GEDI_FILTER <- ONP_CHM_GEDI_FILTER[ONP_CHM_GEDI_FILTER$rh100 > 0,] #lets remove the weird GEDI rh100 values that are 0
mapview(st_as_sf(ONP_CHM_GEDI_FILTER))

Create a regression model between GEDI and ALS CHM

reg <- lm(ONP_CHM_GEDI_FILTER$hoh_2013_dsm_4 ~ (ONP_CHM_GEDI_FILTER$rh100))
summary(reg)
plot((ONP_CHM_GEDI_FILTER$rh100), ONP_CHM_GEDI_FILTER$hoh_2013_dsm_4 , xlab="GEDI rh100", ylab="ALS CHM (ft)", 
     main="Veg Height") 
abline(reg, col="blue")

REQUIRED FIGURES: First screenshot of the filtered GEDI points using mapview and second, a screenshot of the final regression model.

PUTTING IT ALL TOGETHER

In total for the lab so far you should have:

Now for the last part of the lab

ESRM 433:

All your figures and tables need to be fully captioned and submitted in one PDF document. All figures need to be high quality and captions need to describe in brief what the image is of, what data was used, and how it was created.

In addition please explain in 300 words or less:

SEFS 533:

You will include all of the required figures and tables as well as the short essay, same as 433.

In addition, please explain in 300 words or less: